// -----------------------------------------------------------------------------
// M A R I T I M E   R E S E A R C H   I N S T I T U T E   N E T H E R L A N D S
// -----------------------------------------------------------------------------
// Copyright (C) December 14, 2009 MARIN - All rights reserved - http://www.marin.nl
// -----------------------------------------------------------------------------
// Program  : mermaid
// Module   : interactive waves
// File     : MathPulse.cpp
// -----------------------------------------------------------------------------

#include <new>

#ifndef _WINDOWS
#include "mermaid.h"
#endif
#include "FloatUtils.h"
#include "MathGeometry.h"
#include "MathPulse.h"

// -----------------------------------------------------------------------------
// function MathPulse::MathPulse
// -----------------------------------------------------------------------------
/// constructor
/** detailed description */
// -----------------------------------------------------------------------------
MathPulse::MathPulse(const MathGeometry& a_math, double a_dx, double a_dy)
: super(a_dx, a_dy)
, m_geometry(a_math)
{
    // implementation
}


// -----------------------------------------------------------------------------
// function MathPulse::~MathPulse()
// -----------------------------------------------------------------------------
/// destructor
/** detailed description */
// -----------------------------------------------------------------------------
MathPulse::~MathPulse()
{
    // implementation
}


// -----------------------------------------------------------------------------
// function MathPulse::compute
// -----------------------------------------------------------------------------
/// computes pressure pulse for use with Gert Klopman's wave model
// -----------------------------------------------------------------------------
void MathPulse::compute(int *a_px, ///< x offset of resulting pulse array
                        int *a_py, ///< y offset of resulting pulse array
                        Array2D<REAL>* a_pPulse) ///< pulsedepths
{
    // syntactic sugar
    Array2D<REAL>& pulse = *a_pPulse;

    // determine horizontal bounding box of ship
    // simplistic: take longest dimension
    const mml::Vector3d& radii = m_geometry.getRadii();
    REAL r = std::max(radii(0), radii(1));
    mml::Vector3d pos = m_geometry.getPosition();
    double xmin = pos(0) - r - m_dx;
    double ymin = pos(1) - r - m_dy;
    double xmax = pos(0) + r + m_dx;
    double ymax = pos(1) + r + m_dy;

    // determine indices of bounding box corners in waveheight array
    double x, y;
    modf(xmin / m_dx, &x);
    modf(ymin / m_dy, &y);
    int imin = int(x);
    int jmin = int(y);
    modf(xmax / m_dx, &x);
    modf(ymax / m_dy, &y);
    int imax = int(x);
    int jmax = int(y);

    // set dimensions of pulse array to size of bounding box
    pulse.clear();
    int dim1 = imax - imin + 1;
    int dim2 = jmax - jmin + 1;
    pulse.resize(dim1, dim2);

    // determine offsets of pulse array within waveheight array
    *a_px = imin;
    *a_py = jmin;

    // set all depth values in pulse array to 0
    pulse = 0.0;

    REAL pm = m_geometry.getMaximumMeterWaterPressure();

    if (!m_geometry.isEllipse()) // rectangular pulse, no use of angle!
    {
        for (int i1 = 0; i1 < dim1; ++i1)
        {
            for (int i2 = 0; i2 < dim2; ++i2)
            {
                // determine if point is inside pulse rectangle
                REAL xx1 = (i1 + imin) * m_dx; // point in wavegrid to check
                REAL xx2 = (i2 + jmin) * m_dy;
                if (xx1 > pos(0) - radii(0) && xx1 < pos(0) + radii(0) &&
                    xx2 > pos(1) - radii(1) && xx2 < pos(1) + radii(1))
                {
                    pulse[i1][i2] = -pm;
                }
            }
        }
    }
    else // ellipsoid pulse
    {
        // for now, only works for movements and rotations in the horizontal plane
        REAL theta = -m_geometry.getAngle();
        REAL sintheta = sin(theta);
        REAL costheta = cos(theta);
        REAL X1 = pos(0) * costheta + pos(1) * sintheta;
        REAL X2 = -pos(0) * sintheta + pos(1) * costheta;
        for (int i1 = 0; i1 < dim1; ++i1)
        {
            for (int i2 = 0; i2 < dim2; ++i2)
            {
                // determine if point is inside pulse ellipse
                REAL xx1 = (i1 + imin) * m_dx; // point in wavegrid to check
                REAL xx2 = (i2 + jmin) * m_dy;
                REAL XX1 = xx1 * costheta + xx2 * sintheta;
                REAL XX2 = -xx1 * sintheta + xx2 * costheta;
                REAL pr = FloatUtils::square((XX1 - X1) / radii(0)) +
                    FloatUtils::square((XX2 - X2) / radii(1));
                if (pr < 1.0)
                {
                    pr = sqrt(pr);
                    pulse[i1][i2] =
#ifdef TEST
                        -pm * 0.5 * (1.0 + cos(PI * pr));
#else
                        -pm * (1 - FloatUtils::power(-sin(PI * 0.5 * pr), 12));
#endif
                }
            }
        }
    }
}
